#
#
# *------------------------------------------------------------------
# | PROGRAM NAME:  Replication Code for RISP Special Issue
# | DATE: 27 October 2017
# | CREATED BY:  Francesco Visconti
# | PROJECT FILE: The legislative representation of public opinion
# | policy priorities in Italy.
# *----------------------------------------------------------------
# | PURPOSE:               
# | Code to replicate tables.
# | 
# *------------------------------------------------------------------
#

library(foreign)
library(dplyr)
library(psych)
library(plm)
library(tidyr)
library(lmtest)
library(reshape2)
library(sandwich)
library(tseries)

#Clear objects from workspace
rm(list = ls()) # get rid of any existing data 

ls() # view open data sets

#Load data:
dataset <- read.dta("Your/File/Path/Here/ReplicationFiles/dataset.dta")

#Create PLM data object
dataset <- plm.data(dataset, index = c("Issue", "Date"))

####VARIABLES DESCRIPTION#####

#Issue: categorical variable registering the type of issue
#(Crime, Economics, Education, Health Care, Housing, Immigration, Pensions, Taxation, Terrorism, Unemployment)

#Date: factor variable registering the survey date.

#pub_op: continuous variable measuring the aggregated Eurobarometer (EB) Most Important Problem (MIP) 
#answers on the ten issues.

#BillsMPs: continuous variable measuring the legislative agenda of MPs in the semester prior to the survey date.

#BillsMPsMajority: continuous variable measuring legislative agenda of MPs of the majority in the semester prior to the survey date.

#BillsMPsOpposition: continuous variable measuring legislative agenda of MPs of the opposition in the semester prior to the survey date.

#BillsGOVT: continuous variable measuring legislative agenda of government in the semester prior to the survey date.

#Mattarellum: binary variable registering whether the MPs where elected with Mattarellum electoral law or not.

#government: categorical variable registering the government in power.

#ENPS: continuous variable measuring the effective number of parties in parliament.

#ENPC: continuous variable measuring the effective number of parties in government.

#parl_ideol2: continuous variable measuring the average ideology of parliament based on CHES data.

#gov_ideol2: continuous variable measuring the average ideology of government/majority based on CHES data.

##opp_ideol: continuous variable measuring the average ideology of opposition MPs based on CHES data.

####TABLE 1 - Descriptive statistics####

describeBy(dataset[, 3:7], dataset$Issue)

#####TABLE 2 - Effect of public opinion on the legislative agenda of all MPs.#####
#### MPs AGENDA Models #####
#Fixed effects panel using plm and LEVELS with lagged DV#

#Model 1
fixed.levels.lagDV.mps <- plm(BillsMPs ~ lag(BillsMPs) + lag(pub_op) + Mattarellum + parl_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.mps)
#display the fixed effects (constants for each policy area)
summary(fixef(fixed.levels.lagDV.mps))

#Model 2
#MPs agenda with lagged government agenda
fixed.levels.lagDV.mps.govt <- plm(BillsMPs ~ lag(BillsGOVT) + lag(BillsMPs) + lag(pub_op) + Mattarellum + parl_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.mps.govt)

#Arellano Bond Standard Errors
coeftest(fixed.levels.lagDV.mps, vcovHC(fixed.levels.lagDV.mps, method = "arellano"))
coeftest(fixed.levels.lagDV.mps.govt, vcovHC(fixed.levels.lagDV.mps.govt, method = "arellano"))

#PCSE - Beck and Katz
coeftest(fixed.levels.lagDV.mps, vcovBK(fixed.levels.lagDV.mps))
coeftest(fixed.levels.lagDV.mps.govt, vcovBK(fixed.levels.lagDV.mps.govt))

#####TABLE 3 - Effect of public opinion on the legislative agenda of majority MPs.#####
#### MPs MAJORITY AGENDA Models #####
#Fixed effects panel using plm and LEVELS with lagged DV#

#Model 1
fixed.levels.lagDV.majority <- plm(BillsMPsMajority ~ lag(BillsMPsMajority) + lag(pub_op) + Mattarellum + gov_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.majority)
#display the fixed effects (constants for each policy area)
summary(fixef(fixed.levels.lagDV.majority))

#Model 2
fixed.levels.lagDV.majority.opp.govt <- plm(BillsMPsMajority ~ lag(BillsMPsMajority) + lag(BillsMPsOpposition) + lag(BillsGOVT) + lag(pub_op) + Mattarellum + gov_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.majority.opp.govt)

#Model 3
fixed.levels.lagDV.majority.opp <- plm(BillsMPsMajority ~ lag(BillsMPsMajority) + lag(BillsMPsOpposition) + lag(pub_op) + Mattarellum + gov_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.majority.opp)

#Model 4
fixed.levels.lagDV.majority.govt <- plm(BillsMPsMajority ~ lag(BillsMPsMajority) + lag(BillsGOVT) + lag(pub_op) + Mattarellum + gov_ideol2 + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.majority.govt)

#PCSE - Beck and Katz
coeftest(fixed.levels.lagDV.majority, vcovBK(fixed.levels.lagDV.majority))
coeftest(fixed.levels.lagDV.majority.opp, vcovBK(fixed.levels.lagDV.majority.opp))
coeftest(fixed.levels.lagDV.majority.govt, vcovBK(fixed.levels.lagDV.majority.govt))
coeftest(fixed.levels.lagDV.majority.opp.govt, vcovBK(fixed.levels.lagDV.majority.opp.govt))

#Arellano Bond Standard Errors
coeftest(fixed.levels.lagDV.majority, vcovHC(fixed.levels.lagDV.majority, method = "arellano"))
coeftest(fixed.levels.lagDV.majority.opp, vcovHC(fixed.levels.lagDV.majority.opp, method = "arellano"))
coeftest(fixed.levels.lagDV.majority.govt, vcovHC(fixed.levels.lagDV.majority.govt, method = "arellano"))
coeftest(fixed.levels.lagDV.majority.opp.govt, vcovHC(fixed.levels.lagDV.majority.opp.govt, method = "arellano"))


#####TABLE 4 - Effect of public opinion on the legislative agenda of opposition MPs.#####
#### MPs OPPOSITION AGENDA Models #####
#Fixed effects panel using plm and LEVELS with lagged DV#

#Model 1
fixed.levels.lagDV.opposition <- plm(BillsMPsOpposition ~ lag(BillsMPsOpposition) + lag(pub_op) + Mattarellum + opp_ideol + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.opposition)
#display the fixed effects (constants for each policy area)
summary(fixef(fixed.levels.lagDV.opposition))

#Model 2
fixed.levels.lagDV.opposition.maj.govt <- plm(BillsMPsOpposition ~ lag(BillsMPsOpposition) + lag(BillsMPsMajority) + lag(BillsGOVT) + lag(pub_op) + Mattarellum + opp_ideol + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.opposition.maj.govt)

#Model 3
fixed.levels.lagDV.opposition.maj <- plm(BillsMPsOpposition ~ lag(BillsMPsOpposition) + lag(BillsMPsMajority) + lag(pub_op) + Mattarellum + opp_ideol + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.opposition.maj)

#Model 4
fixed.levels.lagDV.opposition.govt <- plm(BillsMPsOpposition ~ lag(BillsMPsOpposition) + lag(BillsGOVT) + lag(pub_op) + Mattarellum + opp_ideol + ENPS, data = dataset, model = "within")
summary(fixed.levels.lagDV.opposition.govt)

#PCSE - Beck and Katz
coeftest(fixed.levels.lagDV.opposition, vcovBK(fixed.levels.lagDV.opposition))
coeftest(fixed.levels.lagDV.opposition.maj, vcovBK(fixed.levels.lagDV.opposition.maj))
coeftest(fixed.levels.lagDV.opposition.govt, vcovBK(fixed.levels.lagDV.opposition.govt))
coeftest(fixed.levels.lagDV.opposition.maj.govt, vcovBK(fixed.levels.lagDV.opposition.maj.govt))

#Arellano Bond Standard Errors
coeftest(fixed.levels.lagDV.opposition, vcovHC(fixed.levels.lagDV.opposition, method = "arellano"))
coeftest(fixed.levels.lagDV.opposition.maj, vcovHC(fixed.levels.lagDV.opposition.maj, method = "arellano"))
coeftest(fixed.levels.lagDV.opposition.govt, vcovHC(fixed.levels.lagDV.opposition.govt, method = "arellano"))
coeftest(fixed.levels.lagDV.opposition.maj.govt, vcovHC(fixed.levels.lagDV.opposition.maj.govt, method = "arellano"))

#####TABLE 5 - Effect of public opinion on the legislative agenda of government.#####
#### GOVERNMENT AGENDA Models #####
#Fixed effects panel using plm and LEVELS with lagged DV#

#Model 1
fixed.levels.lagDV.government <- plm(BillsGOVT ~ lag(BillsGOVT) + lag(pub_op) + gov_ideol2 + ENPC, data = dataset, model = "within")
summary(fixed.levels.lagDV.government)
#display the fixed effects (constants for each policy area)
summary(fixef(fixed.levels.lagDV.government))

#Model 2
fixed.levels.lagDV.government.opp.maj <- plm(BillsGOVT ~ lag(BillsGOVT) + lag(BillsMPsMajority) + lag(BillsMPsOpposition) + lag(pub_op) + gov_ideol2 + ENPC, data = dataset, model = "within")
summary(fixed.levels.lagDV.government.opp.maj)

#Model 3
fixed.levels.lagDV.government.maj <- plm(BillsGOVT ~ lag(BillsGOVT) + lag(BillsMPsMajority) + lag(pub_op) + gov_ideol2 + ENPC, data = dataset, model = "within")
summary(fixed.levels.lagDV.government.maj)

#Model 4
fixed.levels.lagDV.government.opp <- plm(BillsGOVT ~ lag(BillsGOVT) + lag(BillsMPsOpposition) + lag(pub_op) + gov_ideol2 + ENPC, data = dataset, model = "within")
summary(fixed.levels.lagDV.government.opp)

#PCSE - Beck and Katz
coeftest(fixed.levels.lagDV.government, vcovBK(fixed.levels.lagDV.government))
coeftest(fixed.levels.lagDV.government.maj, vcovBK(fixed.levels.lagDV.government.maj))
coeftest(fixed.levels.lagDV.government.opp, vcovBK(fixed.levels.lagDV.government.opp))
coeftest(fixed.levels.lagDV.government.opp.maj, vcovBK(fixed.levels.lagDV.government.opp.maj))

#Arellano Bond Standard Errors
coeftest(fixed.levels.lagDV.government, vcovHC(fixed.levels.lagDV.government, method = "arellano"))
coeftest(fixed.levels.lagDV.government.maj, vcovHC(fixed.levels.lagDV.government.maj, method = "arellano"))
coeftest(fixed.levels.lagDV.government.opp, vcovHC(fixed.levels.lagDV.government.opp, method = "arellano"))
coeftest(fixed.levels.lagDV.government.opp.maj, vcovHC(fixed.levels.lagDV.government.opp.maj, method = "arellano"))

#####TABLE 6 - Effects of legislative agendas on public opinion.#####
#### PUBLIC OPINION AGENDA Models #####
#Fixed effects panel using plm and LEVELS with lagged DV#

#Model 1
#MPs
reg.ldv.mps.6 <- plm(pub_op ~ lag(pub_op) + lag(BillsMPs), data = dataset, model = "within")
summary(reg.ldv.mps.6)
coeftest(reg.ldv.mps.6, vcovHC(reg.ldv.mps.6, method = "arellano"))

#Model 2
#Majority
reg.ldv.majority.6 <- plm(pub_op ~ lag(pub_op) + lag(BillsMPsMajority), data = dataset, model = "within")
summary(reg.ldv.majority.6)
coeftest(reg.ldv.majority.6, vcovHC(reg.ldv.majority.6, method = "arellano"))

#Model 3
#Opposition
reg.ldv.opposition.6 <- plm(pub_op ~ lag(pub_op) + lag(BillsMPsOpposition), data = dataset, model = "within")
summary(reg.ldv.opposition.6)
coeftest(reg.ldv.opposition.6, vcovHC(reg.ldv.opposition.6, method = "arellano"))

#Model 4
#Government
reg.ldv.govt.6 <- plm(pub_op ~ lag(pub_op) + lag(BillsGOVT), data = dataset, model = "within")
summary(reg.ldv.govt.6)
coeftest(reg.ldv.govt.6, vcovHC(reg.ldv.govt.6, method = "arellano"))
